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Abstract. It is a clear unique prediction of the cold dark matter paradigm of cosmological 
structure formation that galaxies form hierarchically and are embedded in massive, extended 
dark halos teeming with self-bound substructure or "subhalos". The amount and spatial 
distribution of subhalos around their host provide unique information and clues on the galaxy 
assembly process and the nature of the dark matter. Here we present results from the Via Lactea 
INCITE simulation, a one billion particle, one million cpu-hour simulation of the formation and 
evolution of a Galactic dark matter halo and its substructure population. 



1. Introduction 

According to the standard model of cosmology our Universe is made up mostly of mysterious 
dark energy (70%) and cold dark matter (CDM, 25%). Ordinary matter (baryons) makes up 
only 5% of the mass density in the Universe today. The dark energy component is homogeneous 
and structure formation is dominated by the CDM component. In the last few years cosmological 
N-body simulations have successfully uncovered how the cold dark matter distribution evolves 
from almost homogeneous initial conditions near the Big Bang into the present highly clustered 
state [H [2]. CDM halos are built up hierarchically in a long series of mergers of smaller halos. 
Early low resolution simulations and simple analytical models found that the end products are 
smooth, triaxial halos. Higher resolution cosmological N-body simulations later modified this 
simplistic picture. The merging of progenitors is not always complete: the cores of accreted halos 
often survive as gravitationally bound subhalos orbiting within a larger host system. CDM halos 
are not smooth, they exhibit a wealth of substructure on all resolved mass scales j3|- Subhalos 
are now believed to host cluster galaxies and the satellite galaxies around the Milky Way. 

In 2007 we received an Innovative and Novel Computational Impact on Theory and 
Experiment (INCITE) grant for one million cpu-hours on Oak Ridge National Laboratory's 
Jaguar Cray XT3 supercomputer, to run a one billion particle simulation of the formation and 
evolution of a Milky Way scale dark matter halo in an expanding cosmological volume. The 
resulting simulation, dubbed Via Lacte(^, has resolved the subhalo population at unprecedented 
detail and closer to the host halo's center than ever before. In the following we describe some 



f Actually this simulation is the second in a series of simulations of ever-increasing resolution and is referred 
to in the literature as Via Lactea IL 



of the technical details that lie behind this computational effort, and present some of the first 
scientific results that have emerged from this study. 



2. Technical Details 

Measurements of the cosmic microwave background by the Wilkinson Microwave Anisotropy 
Probe (WMAP) satellite ^ have given us an exquisite view of the state of the universe a mere 
~ 100,000 years after the big bang. At this time the universe was very simple: a complete 
statistical description of the primordial temperature fluctuations observed by WMAP requires 
only six parameters. This simple statistical description of the early universe allows us to generate 
an initial particle distribution at redshift ~ 100 (only 15 million years after the Big Bang) with 
small linear density perturbations matching the WMAP measurements. Our simulation then 
follows the gravitational collapse of these fluctuations into the highly non-linear clustering regime 
today. 



2.1. The N-body technique 

The dark matter density field is sampled using a finite number N of Lagrangian tracer 
whose evolution is governed by Newton's equations of motion in coordinates co-moving 
expansion of the universe [5] : 
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Here a is the expansion scale factor, 
gravitational potential 
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whose field equation is given by 
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where we are allowing for the contribution of a vacuum energy of density pv = A/(87rG). Using 
the 2"^^^ Friedman equation 
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Because the effective gravitational potential is determined by the density contrast {p — pb), the 
computational error on the gravitational potential is set by how accurately the density field 
is sampled. A larger N will thus result in a more faithful representation of the dark matter 
dynamics, and this has been the driving force leading to ever increasing problem sizes and the 
use of powerful supercomputers for their solution. 



2.2. PKDGRAV, a hierarchical tree code 

PKDGRAV [6], the code that we used for the Via Lactea simulation, is an example of a hierarchical 
tree code, in which the particles are grouped together hierarchically and the gravitational force 
on a particle from each group is approximated by a multipole expansion up to low order. The 
hierarchical grouping is realized in a tree structure, in which the whole computational domain is 
split up into ever smaller sub-regions containing successively smaller numbers of particles, down 
to the leaf nodes of the tree which contain only a small number of particles. The calculation 



of the gravitational force felt by each particle proceeds by "walking" the tree, starting with the 
root node that covers the entire computational domain. A cell-opening criterion is applied to 
each node to determine whether to use its multipole approximation, ending the descent along 
this branch, or if the cells needs to be "opened" and its daughter nodes examined. 

The use of such a hierarchical multipole expansion requires only log N interactions per 
particle, much less than the — 1 separate force calculations required with direct-summation, 
and this allows much larger particle numbers to be used. The accuracy of the force calculation is 
determined by the order to which the multipole expansion is carried out and by how many nodes 
it is applied to. Increasing the order of the expansion allows greater accuracy, but comes at an 
additional computation cost, pkdgrav extends the multipole expansion up to fourth order. 

PKDGRAV employs a spatially balanced binary tree, a relative of the particle balanced k-d tree. 
Starting with the root node, every node is recursively split up into two rectangular pieces of equal 
volume by bisecting the tight bounding box containing all particles in the node. This splitting 
continues down to the leaf nodes, or "buckets", which are allowed to contain at most 16 particles. 
Such a tree has been shown to be asymptotically superior to density based decompositions like 
the k-d tree [7]. 

The cell opening criterion applied by pkdgrav is 

where 6max is the maximum distance from any point in the cell to the cell's center of mass, and 
is a user-specified accuracy parameter typically set to 9 = 0.7 for z < 2 and 9 = 0.5 at higher 
redshifts. If during the walk a cell is deemed sufficiently far from bucket Bi to allow its multipole 
approximation, then it is added to iJj's list of particle-multipole interactions. Alternatively, if 
a particular branch is followed all the way down to its leaf node, then all the particles in that 
bucket are added to -Bj's particle-particle interaction list. 

In order to avoid unphysical two-body collisions the potential is softened for small separations, 
\xi — Xj\ < e. This is accomplished by using a single particle density distribution that is given 
by a Dirac (5-function convolved with a softening kernel, here taken to be a compensating Ki 
kernel At separations greater than about 2e the potential becomes exactly Newtonian. 
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2.3. A new dynamical time step criterion 

An accurate determination of the gravitational potential is a necessary, but not sufficient, 
requirement for a reliable solution to the equations of motion of dark matter. Of equal 
importance is the method used to integrate the equations and the choice of time step. 
Cosmological simulations typically exhibit a dynamic range of densities spanning many orders 
of magnitude. Denser particle distributions require shorter time steps for a given accuracy, but 
using a global time step set by the densest regions in the simulations becomes computationally 
impractical, since the majority of the particles reside in much lower density regions. A more 
efficient use of computational resources can be achieved by the use of individual time steps that 
adapt to the local density. PKDGRAV employs an adaptive leapfrog integrator, which allows such 
adjustable time steps. The time steps are quantized in powers of two of the basic time step Tq. 

A widely used time step criterion is based on a rather ad-hoc combination of the force 
softening and the acceleration a given particle feels, AT < 0.2^J^^. This criterion results in a 

very tight time step distribution and has been shown j9] to produce accurate results at all but 
the very highest densities. However, this criterion is not directly related to the local dynamical 
time, and it has been shown [10] that for very high dynamic ranges (> 10^) it is not quite 
adaptive enough, resulting in overly conservative time steps in the outer low density regions and 
insufficiently short time steps in the very centers of dark matter halos. Instead Zemp et al. [10] 




Figure 1. The distribution of time steps for a density profile of slope 7 = 1, for the standard 
time step AT oc ^ye/\a(r)\ (in red) and for the new dynamical time step AT oc 1/ Gp{< r) (in 
blue). From Zemp et al. (2007) [TO]. 



advocate a time step criterion that is directly set by the local dynamical time a particle feels, 

AT < n , ^ 

y/Gpenc{r) 

where penc is the enclosed density from the particle's location to the dominant dynamical 
structure, and rj is typically set to 0.02 — 0.03. The challenge is to quickly determine the 
dominant structure governing the orbit of a given particle, and Zemp et al. [1^ have achieved 
this by cleverly making use of the hierarchical tree structure used for the force calculation. 
Figure [1] (taken from [10] ) shows a comparison of the time step distribution using the old and 
new time step criteria, and reveals that the new time step is more adaptive, with smaller time 
steps close to the center and larges one in the outer regions. The most recent version of pkdgrav 
contains an implementation of this new algorithm, and it allowed our Via Lactea simulation to 
push to unprecedented central densities at a manageable computational cost. 

An unfortunate side effect of the use of individual time steps is that it destroys the symplectic 
nature of the integrator. This would be highly undesirable for any long term integration of 
planetary orbits, for example. In cosmological simulations, however, dark matter particles 
typically complete a comparatively much smaller number of orbits over a Hubble time, and 
so any secular error is not likely to grow significantly during the simulation. In this case the 
computational advantage of using individual and adaptive time steps outweighs the negative 
side-effects. 
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Table 1. Simulation Parameters 
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40.0 


40.0 


104 1.05 X 10^ 


4.10 X 10^ 


402 


1.93 X 10^2 


201 


60.0 


53,653 



Box size Lbox, (spline) softening length e, initial redshift Zi, number A'hircs and mass Mhircs of high 
resolution dark matter particles, host halo r2oo (radius enclosing 200 times the mean density), M200 
(mass within r2oo), (the peak of the circular velocity curve, Vdrc — \/GM{< r)/r) and rymax (the 

radius at which Vinax occurs), and number of subhalos within r2oo -^sub at z = 0. Force softening lengths 
e are constant in physical units back to z = 9 and constant in co-moving units before. 



2.4- Multi-mass initial conditions 

The Via Lactea simulation makes use of multi-mass initial conditions [111 |9] , which allow the 
halo of interest to be resolved with very high resolution, while at the same time allowing for 
the large scale cosmological tidal field to be represented at much lower resolution. The overall 
computational volume is (40 MpcQ)^, but at the end of the simulation 95% of all high resolution 
particles are found within the inner (2 Mpc)'^. In total the Via Lactea simulation consists of a 
hierarchy of three particle masses: (16,689,261, 28,645,888, 1,048,772,608) particles with masses 
(32,768, 64, 1) times the high resolution particle mass of Mhii-es = 4.10 x 10^ Mo. The remaining 
parameters of the Via Lactea simulation are summarized in Table [TJ 

3. Results 

The main scientific driver for our Via Lactea simulation is to better characterize the dark matter 
subhalo population of a Milky- Way-scale halo. A visual impression of the enormous abundance 
of subhalos is given by Figure EJ which shows a density squared projection of the central (566 
kpc)^ cube. Projecting density squared enhances the contrast between individual subhalos and 
the smooth background. The two insets show progressive zooms into a region of interest. The 
deepest zoom is a projection of a (40 kpc)^ cube centered on subhalo with Vmax = 21.4 km s~B 
and a tidal radius of 46.8 kpc. The high resolution of Via Lactea is able to resolve the second 
level of the subhalo hierarchy - subhalos of a subhalo. 

The overall picture is reminiscent of the "turtles all the way down" image. In fact, in our 
simulation the dumpiness extends over six orders of magnitude in mass down to our subhalo 
resolution limit of ~ 10^ Mq. The theoretical expectation is that this dumpiness continues 
even further, all the way down to the cut-off in the density fluctuation power spectrum set by 
the finite velocity dispersion (or temperature) of the cold dark matter particle [121 113j . some 
10 to 15 orders of magnitude below the scales that we currently resolve. An important caveat 
here is that our simulation does not include the effects of baryons. Interactions with stars or 
giant molecular clouds in the disk of the Milky Way could destroy some fraction of the subhalos 
passing through the center. Such baryonic affects are unlikely to affect the subhalo population 
as whole, however, since most of the subhalos are not on orbits that carry them through the 
central disk. 

3.1. The dark matter subhalo population 

The density profiles of isolated dark matter halos are known to follow a common radial 
dependence, falling off as r~^ in the outer regions, becoming shallower towards the center, and 

^One megaparsec (Mpc) is approximately 3 million light years. 

*Vmax is the peak of the circular velocity curve {vdrc = \/GM{< r)/r) and is used as a proxy for mass. 



Figure 2. A projection of density squared in a 566 kpc x 566 kpc x 566 kpc region centered 
on the host halo at z = 0. The insets (150 kpc and 40 kpc deep, respectively) show zooms into 
a region harboring a subhalo exhibiting abundant sub-substructure. 



finally approaching a central cusp of r~'^, with the exact value of 7 being somewhat uncertain, 
but less than 1.5 [Ml [151 HB]- The very large number of particles in the Via Lactea simulation 
and the extremely accurate integration of their orbits has allowed a precise determination of the 
density profile within the inner kpc of the Galactic halo and its satellite halos, see the left panel 
of Figure [31 The host halo density profile is best fit with a central cusp of slope 7 = 1.24, and 
we find that the inner profiles of subhalos are also consistent with cusps, with slopes scattered 
around 7 = 1.2. At larger radii the subhalo density profiles generally fall off faster than the host 
halo (see inset). The explanation is simply that subhalo density profiles are modified by tidal 
mass loss, which strips material from the outside in, but does not affect the inner cusp structure 




Figure 3. Properties of the subhalo population in Via Lactea, from Diemand et al. (2008) 
[33j . Left panel: Density profiles of the host halo (thick black line) and of eight large subhalos 
(thin colored lines). The profiles remain cuspy down the estimated convergence radius of 380 pc 
(vertical dotted line). The rescaled profiles in the inset demonstrate that i) subhalos are quite 
self-similar in the inner regions and ii) many subhalos are tidally truncated in the outer region. 
Right panel: The number of subhalos above Imax within r2oo = 402 kpc (thick line) and within 
100 and 50 kpc of the halo center (thin lines). Fmax is the peak of the subhalo circular velocity 
fcirc = y/GM{< r)/r, and serves as a proxy for subhalo mass. The distribution is well fit by 
a single power-law, iV(> Kiax) = 0.036(ymax/Kiax,host)~^ (dotted hue). Below Fmax - 3.5 km 
s~^ the subhalo abundance is artificially reduced due to numerical limitations. The inset shows 
the abundance of sub-subhalos in eight large subhalos, and is quite similar to the scaled down 
version of the host halo. 



[nmi]. 

Overall Via Lactea predicts a remarkably self-similar pattern of dark matter clustering 
properties. The right panel of Figure [3] shows that subhalos appear to be scaled down versions 
of their host halo not just in the inner mass distribution, but also in terms of their relative 
substructure abundance. The overall subhalo abundance as a function of mass, as measured by 
the cumulative Knax function, is well described by a single power law 

iV(> Kiax) OC (Fmax/^max,host)"^ (9) 

independently of distance, within 50, 100, or 402 kpc of the center. Furthermore, for the first 
time it has become possible to determine the abundance of subhalos at the second level of the 
substructure hierarchy, and the sub-subhalo Imax-functions (see inset) appear to be quite similar 
to the scaled down version of the host halo's subhalo distribution. 

3.2. The "missing satellites problem" 

It has been known for many years that our host galaxy, the Milky Way, is orbited by a number of 
so-called dwarf galaxy satellites. The original census of about 11 dwarfs [19j has recently doubled 
in size by the discovery of ~ 12 additional ultra-faint dwarfs in the Sloan Digital Sky Survey 
(SDSS) [SniEIlESlEallMlESllMlET]. Many, if not all, of these systems are strongly dark matter 



dominated in their central regions [281I29], and are thus thought to be the luminous counterparts 
to the massive end of the dark matter subhalo population. The apparent discrepancy between the 
number of observed dwarf galaxies (around 50, after correcting for the incomplete sky coverage 
of the SDSS) and the vastly larger number of dark matter subhalos predicted by numerical 
simulations (~ 50,000 in Via Lactea) is well known as the "missing satellites problem". 

The solution to this problem could lie in a modification of the standard cold dark matter 
paradigm of structure formation. For example, warm dark matter, with its higher intrinsic 
velocity dispersion, would lead to a suppression of power at small scales and greatly reduce 
the number of low mass subhalos seen in CDM simulations. On the other hand, the standard 
CDM picture might be fine, if complicated "gastrophysics" , such as radiative or thermal feedback 
from young galaxies, quasars, or supernovae, prevents star formation in the vast majority of dark 
matter subhalos. In this case the Galactic halo would be populated by an enormous number of 
invisible dark matter clumps. The resolution of the "missing satellites problem" depends on both 
progress from the observational side: deeper surveys, more and improved measurements of stellar 
kinematics in the centers of dwarfs, etc.; and from the theoretical side: better characterization of 
the properties of dark matter subhalos from simulations, a better understanding of the efficacy 
of feedback in preventing star formation, etc. Our Via Lactea simulations are one step in this 
direction [30]. 

3.3. Indirect detection of dark matter substructure? 

An intriguing, if somewhat speculative, possibility for learning more about the dark matter 
subhalos of our Milky Way arises if the dark matter particle is a WIMP, a weakly interacting 
massive particle, like the neutralino in the super symmetric extension to the standard model of 
particle physics. In that case the dark matter particle would be its own anti-particle, and in 
regions of sufficiently high density colliding dark matter particles would pair-annihilate. The 
energetic standard model particles that are produced in such annihilations would decay in a 
particle cascade, resulting in a potentially observable signal of electron-positron pairs, neutrinos, 
and gamma-rays. Space based experiments (like PAMELA and the upcoming AMS-02, [3lJ) 
could detect anti-particles produced in dark matter annihilations within about one kpc, and 
huge neutrino telescopes built deep into Antarctic ice (e.g. AMANDA-II and its successor 
IceCube [32j) or in the Mediterranean sea (ANTARES [32]) have also been pursuing such an 
indirect detection. The strength of such a signal depends on the local dark matter density in 
the solar neighborhood, and can be boosted by the dumpiness of the dark matter distribution. 
Extrapolating from the measured subhalo abundance in our simulation, we have estimated a 
local boost (within 1 kpc) of about 40% [33]. 

The hope for the gamma-ray signal rests on ground-based atmospheric Cerenkov telescopes 
(like H.E.S.S., VERITAS, MAGIC, and STACEE, [34J) and on NASA's recently launched 
Gamma-ray Large Area Space Telescope (GLAST, |35]). We have used the Via Lactea simulation 
to make quantitative predictions for the possible dark matter annihilation signal that GLAST 
may detect [36j- By placing fiducial observers at 8 kpc (the distance from the Sun to the Galactic 
center) from the host halo center and summing up the annihilation fluxes from all particles in 
our simulation, we produced allsky maps of the gamma-ray flux from dark matter annihilations 
within the Galaxy. An example of such an allsky map is shown in Figure HI Converting these 
flux maps into expected number counts of gamma-ray photons requires assumptions about the 
nature of the dark matter particle, namely its mass, its annihilation cross section, and the 
gamma-ray spectrum resulting from a single annihilation event. Particle physics theory today 
does not uniquely predict these quantities, and so we have sampled a region of the theoretically 
motivated parameter space. Using detailed models of the known astrophysical gamma-ray 
backgrounds that the dark matter annihilation signal from Galactic substructure has to contend 
with, we determined the expected detection significance for every subhalo in our simulation. 
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Figure 4. Allsky maps (in a Hammer-AitofF projection) of the annihilation signal 
(oc /9^c^ dl) from Via Lactea, for an observer 8 kpc from the halo center. Left panel: the total 
signal from all particles within r2oo- Right panel: The signal from only the subhalo particles. 
From Kuhlen et al. (2008) [36j . 



Although the vast majority of the subhalos are too faint to be detectable, we found that for 
reasonable choices of the particle physics parameters, GLAST should be able to detect a handful 
of individual subhalos at more than five sigma significance. 

4. Summary 

The Via Lactea simulation is a one billion particle simulation of the formation and evolution of 
a Milky Way scale dark matter halo in an expanding cosmological volume. At a cost of about 
1 million cpu hours on ORNL's Cray XT3 supercomputer Jaguar, it provides an unprecedented 
view of the subhalo population in a dark matter halo like the one harboring our own Milky Way. 

PKDGRAV, the code we used for our simulation, is an example of a hierarchical tree code 
that uses a multipole expansion of the gravitational potential to solve the N-body problem in 
a fast and accurate manner. In Via Lactea we employed for the first time a new dynamical 
time step criterion that ensures a highly accurate orbit integration even for particles in very 
dense environments, while at the same time allowing longer time steps for particles in the outer 
regions of the halo. This allowed us to follow the central density profile and resolve the subhalo 
population closer to the host halo center than ever before. 

We find that the dark matter clustering properties are remarkably self-similar: isolated halos 
and subhalos contain about the same relative amount of substructure and both have cuspy inner 
density profiles. The resolution of the "missing satellites problem", a well known discrepancy 
between the number of luminous dwarf galaxies orbiting our Galaxy and the predicted number 
of dark matter subhalos, will require either a modification of the standard cold dark matter 
paradigm of structure formation or baryonic physics that suppresses star formation in all but 
a small fraction of subhalos. In the future it may be possible to indirectly detect dark matter 
substructure, through the measurement of electron-positron pairs, neutrinos, or gamma-rays 
produced in pair-annihilations of dark matter particles in high density regions such as the centers 
of subhalos. 
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